Escaped Muons

Introduction

Part of the Accelerator Readiness Review is looking at muons that leave the ring. To this end, we (Renee, Nathan, James, and myself) have a preliminary version of MDC-2 with full injection, unified fields, and “ghost” detectors in a Geant parallel world (so these sensitive volumes can be placed without overlapping physical structures). We can remove various magnets and fields and see how that affects muons escaping the beam.

# Load libraries
library(readGallery)
library(stringr)
library(dplyr)
library(ggplot2)
library(rgl)
library(readr)
library(purrr)
library(parallel)

The useDataProduct calls should all come together too.

useDataProduct('std::vector<gm2truth::TrackingActionArtRecord>')
useDataProduct('std::vector<gm2truth::RingTrackingPlaneArtRecord>')
useDataProduct('std::vector<gm2truth::GhostDetectorArtRecord>')

Samples

I have run preliminary MDC-2 (using Renee’s FCL files), making 10K events for various scenarios. Files are currently stored in /pnfs/GM2/scratch/users/lyon/arr_20170313 and directories within. Ghost detectors included are Renee’s cylinder just at the outer vacuum wall and encompossing the inflector as well as my ghost detector that is a cylindrical shell just on the inside of the world cube. Note that this code includes the magnet yoke steel.

Find the data

Need a function to turn a /pnfs path into a xrootd url

xrootdify <- function(p) {
  # We /pnfs/BLA --> root://fndca1.fnal.gov/pnfs/fnal.gov/usr/BLA
  paste0('root://fndca1.fnal.gov/pnfs/fnal.gov/usr', str_replace(p, '/pnfs', ''))
}

# Let's look at only the 180* job series (so we just get 10,000 events per scenario)
system("ssh lyon@gm2gpvm04.fnal.gov 'ls /pnfs/GM2/scratch/users/lyon/arr_20170313/180*/*.root'", intern=T) %>% xrootdify() -> arrFiles
arrFiles

[1] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root"   
[2] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089095_0/ARR_unified_noDipole_cyl.root"     
[3] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089097_0/ARR_unified_noInflector_cyl.root"  
[4] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089099_0/ARR_unified_noKickerQuads_cyl.root"
[5] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089123_0/ARR_unified_noKicker_cyl.root"     
[6] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089132_0/ARR_unified_noQuads_cyl.root"      

Let’s pluck out the scenario.

scenarios <- str_match(arrFiles, 'arr_20170313/.+_unified_(.+)_cyl')[,2]
scenarios

[1] "everything"    "noDipole"      "noInflector"   "noKickerQuads" "noKicker"      "noQuads"      

What data are in these file? Well, product_sizes_dumper doesn’t seem to work across XRootD. That’s too bad.

Here’s a picture from Root.

Look at Tracking action data to see how many muons decay and how many escape

Create a reader class (modify outside of this notebook) for TrackingActionArtRecord

readerClassSkel('gm2truth::TrackingActionArtRecord', writeFile = 'trackingActionReader.py')

trackingActionReader

For documentation, here is the TrackingActionArtRecordReader reader class

readr::read_file('trackingActionReader.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class TrackingActionArtRecordReader(GalleryReaderBase):
  def __init__(self, inputTag, evLimit):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID', 
                  'parentTrackID', 'volumeUID', 'status', 'x', 'y', 'z', 'px', 'py', 'pz']
    self.evLimit = evLimit

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))

  def fill(self, ROOT, ev):
    
    if ev.eventEntry() > self.evLimit:
      return False
    
    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2truth::TrackingActionArtRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2truth::TrackingActionArtRecord
      for e in p:                             # Loop over elements and fill
      
        # Let's only accept the muon
        if e.trackID != 1 or e.parentTrackID != 0:
          continue

        self.vals.append(
          [ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, e.parentTrackID,  e.volumeUID, 
            e.status, e.x, e.y, e.z, e.px, e.py, e.pz ])

    return True

Make an instance

trackingActionReaderClass <- createReaderClass_from_file('trackingActionReader.py')$TrackingActionArtRecordReader

Load the data - let’s just look at one file since it takes a long time

taReader <- trackingActionReaderClass(artInputTag("artg4"), 20000)
getGalleryData(arrFiles[1], taReader)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root
Closing file, read 696635144 bytes in 163 transactions
Timings: Overall time = 106.982 s
Time per event: min=0.000   mean=0.009   max=6.802

As of right now, the dCache XRootD door is down - no - it’s back up

taDF <- galleryReader_df(taReader)
nrow(taDF)

[1] 20000

Set the scenario (we’ll need that again, so make it a function)

setScenario <- function(df) {
  df %>% mutate(scenario = factor(fileEntry, levels=0:(length(scenarios)-1), labels=scenarios))
}

taDF %>% setScenario() %>% select(scenario, everything()) -> taDF
taDF

How many track action hits do we get per scenario.

taDF %>% group_by(scenario) %>% tally()

Let’s just get the death of the muons.

Muon deaths in tracking action

taDF %>% filter(status == 1) -> taDeathDF

We need the volume ID. I have a art FCL to run to get this information. I have files and they each may have a different set of volume IDs. Let’s try to do this in parallel. The examples in the help for parallel::mcparallel seem to be useful here.

Create the command strings,

runForVolIDString <- function(i) {
  str_interp(
    "PVS_CSVOUT=${csvout}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
    list( csvout=scenarios[i], fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=arrFiles[i]) 
  )
}
runForVolIDStrings <- sapply(1:length(arrFiles), runForVolIDString)
runForVolIDStrings

[1] "PVS_CSVOUT=everything_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root"      
[2] "PVS_CSVOUT=noDipole_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089095_0/ARR_unified_noDipole_cyl.root"          
[3] "PVS_CSVOUT=noInflector_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089097_0/ARR_unified_noInflector_cyl.root"    
[4] "PVS_CSVOUT=noKickerQuads_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089099_0/ARR_unified_noKickerQuads_cyl.root"
[5] "PVS_CSVOUT=noKicker_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089123_0/ARR_unified_noKicker_cyl.root"          
[6] "PVS_CSVOUT=noQuads_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089132_0/ARR_unified_noQuads_cyl.root"            

Let’s run in parallel

jobs <- lapply(1:length(arrFiles), function(i) mcparallel(system( runForVolIDStrings[i], intern=T ), name=i))
res <- mccollect(jobs)

Now we need to load the csv files. Let’s do that in parallel too

jobs <- lapply(1:length(arrFiles), function(i) paste0(scenarios[i], "_volNames.csv") %>% read_csv(col_names=c("volumeUID", "volName")) %>% mcparallel(name=i) )
volNameTables <- mccollect(jobs)
names(volNameTables) <- scenarios

Let’s see what we got

volNameTables[["everything"]]

Join them together when the scenario name.

volNameTable <- map2_df(volNameTables, scenarios, function(df, sce) df %>% mutate(scenario = sce) )
volNameTable

Now let’s do lookups… [Note that I can’t do what I originally thought, which was split up the data frame into groups by scenario and then put in the factor by group – because in the end the rows get rejoined and then row-binded factor gets messed up.]

taDeathDF -> hold # in case we mess up
taDeathDF %>% inner_join(volNameTable) %>% select(scenario, eventEntry, volName, everything()) -> taDeathDF

Joining, by = c("scenario", "volName", "volumeUID")

taDeathDF

We can check this – we should see different volume names associated with volume IDs for different scenarios

taDeathDF %>% distinct(scenario, volumeUID, volName) %>% arrange(volumeUID)

Let’s plot where things die

taDeathDF %>% group_by(scenario, volName) %>% tally() %>% arrange(scenario, desc(n)) -> deathVolumeTable
deathVolumeTable

Did any of them actually decay?

taDeathDF %>% filter(volName == "ArcSection[00]")

Uh oh!!!! Hardly any decay!

taDeathDF %>% ggplot( aes(x = volName) ) + geom_bar() + theme(axis.text.x = element_text(angle=90, hjust=1))

Look at why hardly any muons decay.

Let’s look at the ‘everything’ scenario.

everythingFile <- which(scenarios == "everything")

I made a new TrackingAction reader

readr::read_file('trackingActionReaderCheckDecay.py') %>% cat

# Store the muon at death and whatever it decayed to at birth

from readGallery import GalleryReaderBase  # Necessary for the base class

class TrackingActionWithDecayReader(GalleryReaderBase):
  def __init__(self, inputTag, evLimit):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID', 
                  'parentTrackID', 'volumeUID', 'status', 'x', 'y', 'z', 'px', 'py', 'pz', 'e', 'turn']
    self.evLimit = evLimit

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))

  def fill(self, ROOT, ev):
    
    if ev.eventEntry() > self.evLimit:
      return False
    
    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2truth::TrackingActionArtRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2truth::TrackingActionArtRecord
      for e in p:                             # Loop over elements and fill
      
        # Let's only accept the muon at death or whatever the muon decays to at birth
        if (e.trackID == 1) or (e.parentTrackID == 1 and e.status == 0):

          self.vals.append(
            [ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, e.parentTrackID,  e.volumeUID, 
              e.status, e.x, e.y, e.z, e.px, e.py, e.pz, e.e, e.turn ])

    return True

Make an instance

taDecayClass <- createReaderClass_from_file('trackingActionReaderCheckDecay.py')$TrackingActionWithDecayReader

Load the data - let’s just look at one file since it takes a long time

taReader <- taDecayClass(artInputTag("artg4"), 20000)
getGalleryData(arrFiles[everythingFile], taReader)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root

tadDF <- galleryReader_df(taReader)
nrow(tadDF)

[1] 52842

Make things easier to look at

tadDF %>% 
  mutate(r    = sqrt(x*x+y*y),
         pmag = sqrt(px*px + py*py + pz*pz) ) %>% 
  select(trackType, volumeUID, status, r, e, pmag, trackID, turn, x, y, z) -> tadDF

So I see a positron in alomst every case. But the muon appears to stop in the iron. Let’s just plot x,z

taDeathDF %>% ggplot(aes(x=z, y=x)) + 
  geom_point() + xlim(c(-10000, 10000))

tmp <- open3d()

ring <- cylinder3d( center=rbind(c(0,-90, 0), c(0,90,0)),
                    radius=7112,
                    sides=20, closed = F)

clear3d()
with(taDeathDF, 
     plot3d(x=x, y=y, z=z, type='p', main='Muon death positions', 
       xlab="x", ylab="y", zlab="z", xlim=c(-9000, 9000), zlim=c(-9000, 9000)))
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

This really doesn’t look like muons are being stored.

What fraction of muons have zero momentum when they die (so they die in something dense),

tadDF %>% 
  filter( trackID == 1, status == 1) %>% 
  summarize(pZeroPct = sum(pmag<0.1)/nrow(.)*100) -> pZeroPct
pZeroPct

So % of muons die in something dense. How many decay?

tadDF %>% filter(status == 1, trackID == 1, pmag>0.1)  -> tad.muonsWithP.DF
tad.muonsWithP.DF

clear3d()
with(tad.muonsWithP.DF, 
     plot3d(x=x, y=y, z=z, type='p', main='Muon death with momentum positions', 
       xlab="x", ylab="y", zlab="z", xlim=c(-9000, 9000), zlim=c(-9000, 9000)))
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

These really don’t look stored.

Ring tracking planes

Look at the RingTrackingPlaneArtRecord.

readerClassSkel(extractClass = "gm2truth::RingTrackingPlaneArtRecord", fillClass = "gm2truth::BasicArtRecord", 
                writeFile = 'ringTrackingPlaneReader.py')

readr::read_file('ringTrackingPlaneReader.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class RingTrackingPlaneArtRecordReader(GalleryReaderBase):
  def __init__(self, inputTag):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'particle',
                  'trkid', 'parentid', 'pdgid', 'nturn', 'fracturn', 'detname', 
                  'detnum', 'time', 'energy', 'x', 'y', 'z', 'px', 'py', 'pz',
                  'radius', 'angle']
    # !!!! Modify the names accordingly - MUST MATCH self.vals.append CALL BELOW

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.RingTrackingPlaneArtRecord))

  def fill(self, ROOT, ev):

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2truth::RingTrackingPlaneArtRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2truth::BasicArtRecord
      for e in p:                             # Loop over elements and fill
        
        f = e.basicInfo

        self.vals.append(
          [ ev.fileEntry(), ev.eventEntry(), f.particle, f.trkid, f.parentid, f.pdgid, f.nturn, e.fturn, f.detname,
            f.detnum, f.time, f.energy, f.pos.x(), f.pos.y(), f.pos.z(), 
            f.mom.x(), f.mom.y(), f.mom.z(), f.radius, f.angle ])

    return True

rtpC <- createReaderClass_from_file('ringTrackingPlaneReader.py')$RingTrackingPlaneArtRecordReader

Make one and load the data

rtpR <- rtpC(artInputTag('artg4:RingTrackingPlanes'))

getGalleryData(arrFiles[everythingFile], rtpR) -> timings

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170313/18089092_0/ARR_unified_everything_cyl.root

How many events did we see

length(timings$eventTimes)

[1] 10000

rtpDF <- galleryReader_df(rtpR)
nrow(rtpDF)

[1] 4

Only 4 events had a ring tracking plane hit?

rtpDF

Nathan fixes it

Nathan pushed some now code that supposedly fixes this. Let’s look.

Ring tracking planes

myNewFile <- '/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root'

getGalleryData(myNewFile, rtpR) -> timings

Opening first file...
Successfully opened file /home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root
Closing file, read 765054 bytes in 57 transactions
Timings: Overall time = 0.130 s
Time per event: min=0.000   mean=0.000   max=0.016

length(timings$eventTimes)

[1] 1000

newRtpDF <- galleryReader_df(rtpR)
nrow(newRtpDF)

[1] 3203

newRtpDF

Now I see a lot more muons going around the ring!!

How many turns?

newRtpDF %>% ggplot(aes(x=nturn)) + geom_histogram(bins=100)

Not sure if this is right, but it looks better than before!

Let’s look at the volume names and see how mnay of these muons decay. It should be around 10% or so.

Tracking action and volume ids

Let’s load the tracking action data and the volume name.

First, load the tracking action data…

tanReader <- trackingActionReaderClass(artInputTag("artg4"), 20000)
getGalleryData(myNewFile, tanReader)

Opening first file...
Successfully opened file /home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root
Closing file, read Timings: Overall time = 2.840 s
66328385 bytes in 597 transactions
Time per event: min=0.000   mean=0.003   max=0.685

tanDF <- galleryReader_df(tanReader)
tanDF

Just look at the muon death

tanDF %>% filter(trackID == 1, status == 1) -> tanDeathDF
tanDeathDF

Put in volume name

runVolString <- str_interp("PVS_CSVOUT=${csvout}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
    list( csvout='nathanFixed', fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=myNewFile) )
runVolString

[1] "PVS_CSVOUT=nathanFixed_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 /home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"

system(runVolString)

%MSG-i MF_INIT_OK:  gm2 21-Mar-2017 14:31:10 UTC JobSetup
Messagelogger initialization complete.
%MSG
21-Mar-2017 14:31:11 UTC  Initiating request to open input file "/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"
21-Mar-2017 14:31:12 UTC  Opened input file "/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"
Begin processing the 1st record. run: 1 subRun: 0 event: 1 at 21-Mar-2017 14:31:12 UTC
21-Mar-2017 14:31:12 UTC  Closed input file "/home/me/gm2/renee_arr/try/mdc2-1/ARR_unified_everything_cyl.root"

TrigReport ---------- Event  Summary ------------
TrigReport Events total = 1 passed = 1 failed = 0

TrigReport ------ Modules in End-Path: end_path ------------
TrigReport  Trig Bit#    Visited     Passed     Failed      Error Name
TrigReport     0    0          1          1          0          0 pvsToFile

TimeReport ---------- Time  Summary ---[sec]----
TimeReport CPU = 0.000000 Real = 0.000137

Art has completed and will exit with status 0.

volNames <- read_csv('nathanFixed_volNames.csv', col_names=c("volumeUID", "volName"))

Parsed with column specification:
cols(
  volumeUID = col_integer(),
  volName = col_character()
)

volNames

Let’s join!

tanDeathDF %>% inner_join(volNames) %>% select(eventEntry, volName, everything()) -> tanDeathDF

Joining, by = "volumeUID"

tanDeathDF

tanDeathDF %>% distinct(volumeUID, volName) %>% arrange(volumeUID)

tanDeathDF %>% group_by(volName) %>% tally() %>% arrange(desc(n))

So, about 1% decay!

tanDeathDF %>% 
  mutate(r    = sqrt(x*x+z*z),
         pmag = sqrt(px*px + py*py + pz*pz) ) %>% 
  select(eventEntry, trackType, volName, status, r, pmag, trackID, x, y, z) -> tanDeathdDF

tanDeathdDF %>% filter(volName == 'ArcSection[00]')

Do these correspond to the muons that make many turns? We need to do a join with the ring tracking planes

tanDeathdDF %>% inner_join(newRtpDF, by="eventEntry", suffix=c(".trk", "rtp")) -> tanDeathRTPDF
tanDeathRTPDF

Let’s look at guys that have more than one turn. Remember, we get a row for every hit. If we really want to see how many turns we get, we need to look at the last hit for the event.

tanDeathRTPDF %>% filter(fracturn >= 1) %>% arrange(eventEntry, fracturn) %>% group_by(eventEntry) %>% do(tail(., 1)) %>%
  select(eventEntry, volName, r, pmag, fracturn )

Hmm - we’re missing some muons. Let’s look at those.

tanDeathRTPDF %>% filter(volName == 'ArcSection[00]') %>% arrange(eventEntry, fracturn) %>% group_by(eventEntry) %>% do(tail(., 1)) %>%
  select(eventEntry, volName, r, pmag, fracturn )

What happened to event 791?

newRtpDF %>% filter(eventEntry == 791)

Let’s look at some muon tracks. How about 791 and 893. Since our file is on disk, this will be faster.